import numpy as np;
import numpy.linalg as lin;
import matplotlib.pyplot as mp;


#L'equation liant H, U, et  V est : 

#H(U).V = -F(U)


#Xk+1 = Xk - J(Xk)-1 F(Xk)

######### Part 1 : Newton-Raphson method

# Methode de Newton Raphson pour trouver racines d'une fonction de dimension 1
def Newton_Raphson(f, df, x0, eps, Nmax):
	j=0;
	x=x0;
	while (j<Nmax and np.abs(f(x))>eps):
		x = x - f(x)/df(x);
		print x;
		j = j+1;
	return x;

# Methode de Newton Raphson avec backtracking	
def Newton_Raphson_B(f, df, x0, eps, Nmax):
	j=0;
	x=x0;
	while (j<Nmax and np.abs(f(x))>eps):
		nf = f(x);
		step = - f(x) / df(x);
		while np.abs(f(x + step)) > np.abs(f(x)):  # Backtracking
				step = step * 2 / 3; 
		x = x + step;
		print x;
		j = j+1;
	return x;
	
# Methode de Newton Raphson pour trouver racines d'une fonction de dimension N
# Avec Backtracking
# F fonction a traiter
# J Matrice Jacobienne de F
# U0 Vecteur initial pour la recherche
# eps ecart minimal a zero recherche
# Nmax nombre maximum d'iteration de l'algorithme
def Newton_Raphson_N(F, J, U0, eps, Nmax):
	j=0;
	U=U0;
	
	while (j<Nmax and lin.norm(F(U))>eps):
		nf = lin.norm(F(U));
		#print U;
		V = lin.lstsq(J(U), -F(U))[0];
		while lin.norm(F(U + V)) > lin.norm(F(U)):  # Backtracking
			V = V * 2 / 3; 
		U = U + V;
		#print U;
		j = j+1;
	return U;
	
###################
# Partie de test : 

fonction1  = lambda x: np.cos(x) - x*x*x;
fonction2  = lambda x: x*x*x;
dfonction1 = lambda x: -np.sin(x) - 3*x*x;
dfonction2 = lambda x: 3*x*x;
	
# fonction de n=4 dimension
def fn(u):
	n=4;
	resultat = np.zeros(n);
	resultat[0] = fonction1(u[0]);
	resultat[1] = fonction1(u[1]);
	resultat[2] = fonction1(u[2]);
	resultat[3] = fonction1(u[3]);
	return resultat;
	
# Jacobien de fn
def dfn(u):
	#On n'a rempli que la diagonale de la matrice !!!
	n=4;
	resultat = np.zeros((n,n));
        #Regarder si les autres derivees sont nulles.
	#Les autres derivees sont nulles car pour tout (i,j), dFi/dxj = 0 ou F est la fonction plus haut en parametre pour le test.
	resultat[0, 0] = dfonction1(u[0]);
	resultat[1, 1] = dfonction1(u[1]);
	resultat[2, 2] = dfonction1(u[2]);
	resultat[3, 3] = dfonction1(u[3]);
	return resultat;

print "################################"
print "# Partie 1: Methode de Newton-Raphson"
print "#"
print "################################"
print ""
print "Test de la methode de Raphson avec la fonction de dimension 1 cos(x) - x^3";
print "Calcul racine avec pour valeur initiale '-1', eps=0.001, 50 iteration max:";
print Newton_Raphson(fonction1, dfonction1, -1, 0.001, 50);
print ""

print "Calcul identique avec Backtracking :";
print Newton_Raphson_B(fonction1, dfonction1, -1, 0.001, 50);
print ""

print "################################"
print "Test de la methode de Raphson avec une fonction de dimension n=4, fi(x)=cos(xi) - xi^3";
print "Calcul racine avec pour vecteur initial [-1,-1,-1,-1], eps=0.001, 50 iteration max:";
print Newton_Raphson_N(fn, dfn, np.array([-1, -1, -1, -1]), 0.001, 50);

print "Calcul racine avec pour vecteur initial [-1,-2, 3, 1], eps=0.001, 50 iterations max:";
print Newton_Raphson_N(fn, dfn, np.array([-1, -2, 3, 1]), 0.001, 50);

#Remarque : on pourrait faire des graphes qui montrent la vitesse de convergence selon eps et le nombre d'iterations.
#Et selon le vecteur d'origine !!! (pas pareil [-1,-1,-1,-1] que [-1, -2, 3, 1])

def graph_tests_nb_iterations(f, df, eps):
	s = np.arange(1,11,1)
	t = np.arange(0,10,1)
	
	for i in np.arange(0,s.size):
		t[i] = lin.norm(Newton_Raphson_N(f, df, np.array([-1,-2, 3, 1]), eps, i));
		print t[i]
	mp.clf();

	g = mp.plot(s, t, linewidth=1.0);

	mp.xlabel('Nombre iterations de la methode Newton-Raphson')
	mp.ylabel('Norme du vecteur')

	mp.savefig('Convergence de la methode Newton-Raphson.png')

graph_tests_nb_iterations(fn, dfn, 0.001)
graph_tests_nb_iterations(fn, dfn, 1)
